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ABSTRACT 


The  estimation  of  the  Doppler  shift,  frequency  scale  factor  and 
amplitude  of  the  spectrum  of  a  stationary  Gaussian  random  process 
from  a  record  of  limited  duration  is  considered,  ’  ’ith  the  approximation 
that  the  coefficients  of  the  Fourier  series  expansion  of  a  realization 
of  long  time  duration  are  independent,  maximum  likelihood  estimates  are 
derived  and  the  Cramer -Rao  lower  bounds  to  the  variances  are 
evaluated. 
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I. 


Introduction 


The  problem  of  estimating  parameters  of  the  spectrum  of  a  stationary  random 
process  from  a  record  of  limited  duration  arises  in  a  number  of  applications.  For 
example,  in  radar  astronomy  the  signal  returned  from  the  moon  or  a  planet  is  a  random 
process.  The  radial  velocity  of  such  a  body  can  be  inferred  from  the  Doppler  shift  and 

T 1  2l 

the  rotation  rate  from  the  spectral  width. L  ’  J  Considerable  attention  has  been  directed 
toward  estimation  of  over-all  spectral  shapes.^  ^  However,  for  situations  like  the 
above  in  which  the  spectral  shape  is  assumed  to  be  known  except  for  one  or  more 
parameters  the  treatment  has  been  quite  limited.  The  only  reference  we  can  cite  is  that 
by  Swerling^  who  has  determined  an  estimate  of  the  frequency  scale  factor  h  which 
is  optimum  among  a  certain  class.  He  first  forms  an  over-all  spectral  estimate  <fc*(f) 

r  41 

from  the  data  by  the  methods  of  Blackman  and  Tukey  J  and  then,  assuming  a  knowledge 
of  the  approximate  value  of  h  and  Gaussian  statistics,  selects  as  an  estimator  that  linear 
functional  of  4>*(f)  subject  to  certain  restrictions  which  has  a  minimum  mean  square 
error.  His  expression  for  the  mean  square  error  is  found  to  be  approximately  equivalent 
to  (21)  or  (27)  below  when  the  spectrum  is  symmetrical  about  some  center  frequency  and 
<fr*(f)  is  estimated  at  a  large  number  of  points,  and  is  generally  larger  otherwise. 

II.  The  Estimation  of  Doppler  Shift 

r  7  s] 

The  approach  followed  here  is  to  apply  the  method  of  maximum  likelihood1  ’  J 
with  certain  approximations.  The  following  assumptions  are  made: 

a)  A  realization  of  the  random  process  x(t)  is  observed  for  0  s  t  s  T. 

b)  The  process  is  stationary,  Gaussian,  and  zero -mean. 

c)  The  single -sided  spectrum  of  the  process  is 

*(f)  =  P(f-fd) 
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where  P(f-f^)  is  known  except  for  the  translation  parameter  f^,  the 
Doppler  shift,  which  is  to  be  estimated.  If  the  process  consists  of 
a  signal  plus  an  independent  noise,  $(f)  is  merely  the  sum  of  the 
two  individual  spectra. 

d)  f ,  lies  within  a  known  finite  interval,  f  ,  .  —  f ,  —  f. 

d  dmin  d  dmax 

e)  For  any  value  of  f^  within  this  interval  only  values  of  4>(f)  within  a 

finite  range  0  <  f ^  <  f  <  f^  are  dependent  upon  fj.  f^  and  f^  are 

independent  of  f  and,  for  convenience,  are  integral  multiples  of  1/T. 

This  assumption  eliminates  unimportant  end  effects  from  the 

development  and  is  valid  if  P(f-f  )  has  some  constant  value  $  out- 
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side  of  a  fixed  frequency  range  or  if  x(t)  is  band-limited  after  the 
Doppler  shift  has  been  imparted.  See  Figure  1. 

f)  P(f-f^)  obeys  certain  general  regularity  conditions  so  that  it  is  free 
of  discontinuities,  varies  slowly  over  any  interval  of  width  1/T  and  is 
non-zero  for  f  sfs  f^. 

The  likelihood  function  of  the  realization  is  now  expressed  in  terms  of  the  Fourier 
coefficients 
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where 


f  =  1/T 
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(4) 


2 


These  coefficients  are  taken  as  the  "observable  coordinate s’,L  J  of  the  process.  This 

is  plausible  since  under  general  conditions  x(t)  can  be  represented  for  almost  every 
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sample  function  by  the  limit  in  the  mean  of  its  Fourier  series  expansion.  J  The  a^ 
and  b^  have  a  joint  Gaussian  distribution^-  ^  and  with  the  above  definitions 
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For  large  T  these  coefficients  approach  independence.  See  [9] ,  especially  Theorem  7, 
The  principal  approximation  introduced  here,  upon  which  all  following  equations  depend, 
is  that  they  are  independent  for  finite  T.  This  is  quite  accurate  when  assumption  f)  holds. 
Then  with  assumption  e)  from  the  multivariate  Gaussian  distribution  the  logarithm  of 
the  likelihood  function  becomes 


* '  1 

where  =  Tf^,  =  Tf^,  and  C  with  a  subscript  denotes  a  quantity  independent  of  the 
parameter  to  be  estimated.  The  first  summation  is  independent  of  f  to  an  excellent 
approximation  from  assumption  f).  Hence  the  maximum  likelihood  estimate  f  can  be 
approximated  by  minimizing 
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+  C, 


(6) 


N, 


A(fH)  =  E 


M  P(nf  *  O 

n=N^  0  d 


(7) 


This  can  be  thought  of  as  a  cross -correlation  between  the  observed  spectrum  c^  and  the 
reciprocal  of  P(nfQ-  f^).  Knowledge  of  the  absolute  amplitude  level  of  the  spectrum  is 
not  required. 
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A  necessary  condition  that  a  maximum  likelihood  estimate  must  satisfy  is 
£’(fd)  =  0  where  the  prime  denotes  differentiation  with  respect  to  the  argument.  Applied 
to  (7),  this  shows  that  f  must  satisfy 

^2  P'(nf  ■  f  ) 

E  c l  -2-°-d-  =0  (8) 

n=Ni  P  <"V  fd> 

When  the  approximate  minimum  of  (7)  is  found  the  determination  of  the  point  at  which 
this  cross -correlation  passes  through  zero  locates  the  minimum  exactly. 


P(f"fd)  =  4>o  +  p(f'fd)  (9) 

where  max  J  p(f  -  f^)  |  «  4>q  then 

'A<fd>  =  ~2  S  \  P<"V'd>+C2  (10) 

$  n 
o 

and  f  is  given  approximately  by  maximizing  this  summation. 
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In  calculating  f  nearby  values  of  can  be  averaged  providing  that  $(f)  does  not 
vary  appreciably  over  the  included  interval.  Therefore  A(f^)  can  be  approximated  by 

T/2-^t>d(  <»> 

fi  d 

where  «fc*(f)  is  a  smoothed  spectral  estimate  obtained  by  the  Blackman -Tukey  method  or 

2 

by  a  spectrum  analyzer  as  long  as  all  of  the  c^  are  effectively  included  and  the  detailed 
structure  of  the  spectrum  is  not  obscured. 

The  sampling  errors  of  these  estimates  are  now  investigated  by  applying  the 


following  general  results  which  hold  under  suitable  regularity  conditions.  For  a  parameter 


5 


since 


B  cn  *  2fK  '  fd> 

and 

E  c4  =  8P2(nf  -  f .) 
n  o  d 

With  fj  set  equal  to  its  true  value 


E[£'(fH)] 


4>'(nf  ) 
_ 0_ 

^(nf  ) 
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and  for  the  entire  set  of  coefficients 
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(20) 


Therefore,  from  the  asymptotic  properties  of  maximum  likelihood  estimates,  for  large 

T,f ,  tends  to  a  Gaussian  distribution  with  mean  f ,  and  variance 
d  d 


(21) 


where  the  summation  (20)  has  been  approximated  by  an  integral.  The  variations  of  the 
structure  of  the  spectrum  which  permit  estimation  of  f^  are  thus  measured  by  the 
integral  in  (21).  The  expression  (21)  is  also  the  Cramer -Rao  lower  bound  to  the  mean 
square  error  of  any  estimate,  providing  db/df^  =  0,  which  will  be  true  of  any  spectrum 
cross -correlation  type  of  estimate  from  symmetry  considerations. 

As  a  digression,  an  interesting  phenomenon  concerning  (21)  is  pointed  out.  In 
general,  the  Cramer-Rao  lower  oound  (12)  for  the  variance  of  an  unbiased  estimate  can 
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be  written  as 


2  V1  f  -  f'te*<0  2 
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where  x  is  the  vector  random  variable  of  the  observations  and  f(x)  is  the  probability 
density.  This  is  almost  identical  in  form  to  (21)  although  f(x)  and  <fc(f)  represent  entirely 
different  quantities.  This  similarity,  which  has  also  been  found  to  hold  for  other 
functions  of  £  in  this  situation,  is  entirely  unexpected  and  does  not  occur  in  other 
parameter  estimation  problems. 

III.  Estimation  of  Frequency  Scale  Factor  (Spectral  Width) 

The  assumptions  of  Section  II  are  modified  as  follows: 
c’)  The  single  sided  spectrum  of  the  process  is 

*(f)  =  P(h(f-f  )) 

c 

where  f^  is  a  known  center  frequency  and  P(h(f  -  f^))  is  known 

except  for  the  frequency  scale  factor  h  which  is  to  be  estimated. 

d’)  h  lies  within  a  known  finite  interval  0  <  h  ^  h  ^  h 

min  max 

e')  and  f)  are  the  same  as  e)  and  f)  except  that  f  is  replaced  by  h. 

Sc:e  Figure  2. 

Note  that  if  f  =  0  and  the  non-constant  part  of  <fc(f)  is  confined  to  a  narrow  high- 
frequency  band  then  the  estimation  of  h  is  approximately  equivalent  to  that  of  f  . 


Proceeding  as  before ,  h  is  obtained  by  maximizing 


£<h)  =  -  ~  Yj  2l°S  1*  P(h  <nfc  -  fc))  +  Y 


P(h(nf  -f  ))  +C4 
n=Nj  o  c 


/  2log2.P(h(f-f  »df+  /  ..  df  +  C 


'  P(h  (f  -  f  » 
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The  evaluation  of  Ij  can  be  simplified  by  writing 

P(h(f-fc))  =  *o[l  +p(h(f-fc))] 


Then 


I  =  2(f  -f,)  2  log  2tt  <t>  +/  2  log  2n  [  1  +  p(h  (f  -  f  ))]  df 

12  1  o  f  '' 

1 

=  C  +T  (24) 
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On  the  basis  of  the  previous  assumptions  the  entire  non-zero  part  of  the  integrand  of 
If  is  always  included  within  the  limits  of  integration  (see  Fig.  3).  By  a  simple 


transformation  of  variables 


ri=h  /  2  log  2tt  [  1  -4-  p(f -fc)] 


Hence 


m  *  -  ^7 


♦*(0  ^  +c 
P(h(f-fc))  df  +C6 


The  asymptotic  properties  of  maximum  likelihood  estimates  previously  mentioned 
show  that  as  T  becomes  large  h  tends  to  a  Gaussian  distribution  with  mean  h  and  variance 


T  /  (f-f  )2  [*’(f)/*(f)]2  df 


which  is  calculated  in  the  same  way  as  (21).  The  effect  of  spectrum  irregularities  is 
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thus  emphasized  by  the  quadratic  factor  (f  -  f£)  .  For  estimates  of  h,  db/dh  may  not 
be  zero  so  the  Cramer -Rao  lower  bound  is  not  so  informative  as  in  the  previous  case. 
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IV.  Estimation  of  Spectrum  Amplitude^ 

It  is  interesting  to  investigate  the  estimation  of  spectrum  amplitude  by  the  same 
method.  Assumptions  a)  and  b)  of  Section  II  are  retained  and  the  following  are  adcJed: 
c”)  4>(f)  =  AP(f)  +  N 


d") 


where  P(f)  and  N  are  known  and  A  is  to  be  estimated, 
o 

The  Fourier  coefficients  of  x(t)  can  be  determined  only  over  the 


finite  frequency  range  f  ^  f  s  f  . 


Proceeding  as  before,  there  is  obtained 

N, 
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and 
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£'<A)  =  -  |  S 
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P(nf  ) 
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AP(nf  )  +  N 
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2[AP(nfo)  +  NQJ 


(29) 


For  any  unbiased  estimate  A*,  the  Cramer-Rao  lower  bound  is  found  to  be 
VarA*  1 


N2 
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(30) 


or  approximately 


Var  A* 


(31) 
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AP(f) 


AP(f)  +  N 


df 


^This  problem  was  brought  to  the  author's  attention  by  Prof.  W.  L.  Root. 
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A  general  explicit  solution  for  the  maximum  likelihood  estimate  A  is  not  known. 
However  there  are  two  significant  special  cases.  The  first  is  for  N  =0.  Then 
setting  £’(A)  =  0,  there  is  obtained 


A  = 


T«2  -  V 


r 2 

Jf  P(f) 


df 


(32) 


This  is  an  unbiased  estimate  and  the  variance  is  found  directly  to  be 

Var  A  1 


1  <f2  "  V 


(33) 


Therefore  as  (f  -  f^)  approaches  infinity  the  variance  approaches  zero  for  any  fixed  T. 
This  result  would  be  expected  from  the  well  known  singularity  of  this  situation. 

A  second  special  case  occurs  when  AP(f)  «  for  all  f.  A  can  be  approximated 
by  setting 


so  (29)  becomes 


AP(nf  )  +  N  «  N 
o  o  o 


£'(A)*  T 


P(nfo) 

N2 


AP(nf  )  +  N 
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(34) 


and 


f  [**(f)  -  Nj  P(f)df 

A  =  -  (35) 

/  P2(f)  df 
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The  estimate  is  also  unbiased  and  gives 


Var  A  _ 


N 


(36) 


T  /  [  AP(f)]  df 
f. 


which  agrees  with  known  radiometer  formulas.  The  variances  (33)  and  (36)  are  found 
to  be  the  same  as  the  Cramer-Rao  lower  bounds  so  these  estimates  are  efficient. 

V.  Application  to  Radar  Astronomy  Measurements 

The  preceding  analysis  applies  directly  to  the  situation  where  a  target  is 
illuminated  by  a  long  pulse  and  the  parameters  of  the  spectrum  of  the  return  signal  are 
to  be  estimated.  It  also  applies  to  a  dual  situation  in  which  the  range  or  extent  of  a 
target  is  to  be  estimated.  Consider  a  target  with  a  time  delay  tq  such  that  when  it  is 
illuminated  by  a  unit  impulse  the  reflected  signal  r(t  -  r  )  is  a  Gaussian  random  process 
with  mean  zero  and 

E  r(t  -  tq)  r(tj  -  tq)  =  R(t  -  tq)  6(t  - 1^ 

so  that  the  received  signal  is  a  non -stationary  white  noise  process.  This  model  is 
0  consistent  with  the  representation  of  the  target  as  a  collection  of  a  large  number  of 
point  scatterers  randomly  and  independently  distributed  with  a  spatial  density  cor- 
responding  to  *7 R(t  -  tq)  . 

Suppose  that  R(t  -  t^)  is  known  except  for  the  translation  parameter  tq  which  is 
to  be  estimated.  The  following  additional  assumptions  are  made: 

a)  The  target  is  illuminated  by  a  unit  impulse  bandlimited  to 

O^f  <  W.  This  assumption  is  necessary  to  represent  the  practical 
situation  and  to  avoid  a  degenerate  result.  If  an  RF  pulse  has  been 
transmitted  the  received  signal  is  assumed  to  have  been  translated 
to  zero  frequency. 
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b)  To  lies  within  a  known  infinite  interval. 

c)  R(t  -  r  )  is  of  finite  duration  and  lies  within  the  interval  t,  stsr. 

o  12 

for  any  allowed  value  of  r  .  For  convenience  take  t,  =  N./2W  and 

o  11 

=  N^/2W  where  and  are  integers. 

d)  R(t  -  t  )  is  slowly  varying  over  any  interval  At  =  1/2W  and  has  no 
discontinuities.  (Unfortunately  this  assumption  is  not  fulfilled  for 
many  physical  targets.) 

e)  The  return  signal  is  perturbed  by  additive  stationary  white  Gaussian 
noise  of  single -sided  spectral  density  Nq. 

f)  The  return  signal  plus  noise  are  passed  through  an  ideal  bandpass 
filter  0  ^  {  s  W  to  give  a  received  signal  s(t). 

On  the  basis  of  these  assumptions  the  bandlimited  signal  s(t)  can  be  represented 
by  2(t  -t  )  W  +  1  sampled  values  s(nAt).  To  a  good  approximation  the  s(nAt)  are 
independent  and 

Var  s(nAt)  =  [R(nAt  -  t  )  +  Nq  ]  W  (37) 

The  logarithm  of  the  likelihood  of  the  set  of  samples  is,  with  this  approximation. 


2  2 
J.  y  _ s  (nAt) 

2  ",  [  R(nAtJ  -  t  )  +  N  ]  W 

n=N,  1  '  o  oJ 


+  C, 


(38) 


The  maximum  likelihood  estimate 
Therefore  as  W  becomes  large  the 


t  is  seen  to  be  analogous  to  the  estimate  f , 
o  d 

variance  of  the  asymptotic  distribution  of  tq 


of  Section  n. 
is 


R*(t  -  T  ) 
_ Q 

R(t  -  T  )  +  N 
o  o 


(39) 
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If  the  effective  duration  of  the  band -limited  impulse  is  taken  as  T  =  1/W  the  variance 
becomes 
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f- 

~i 

0 

(40) 
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f- 
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R(t  -  r  )  +  N 

dt 

I 

l 

L  0  °J 
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Since  the  exact  transmitted  pulse  shape  should  not  greatly  influence  the  final  result, 
this  expression  could  be  taken  as  an  order  of  magnitude  approximation  for  other 
pulse  shapes. 

The  depth  (spatial  extent)  of  a  distributed  target  is  a  scale  parameter  and  there¬ 
fore  can  be  estimated  in  a  manner  exactly  analogous  to  that  of  Section  III. 

VI.  Conclusions 

The  approximations  which  have  been  employed  are  of  increasing  accuracy  as 
the  observation  time  becomes  large  and  seem  to  lead  to  useful  results.  However,  a 
further  quantitative  justification  would  certainly  be  desirable.  The  joint  estimation  of 
several  parameters  is  a  straightforward  extension  of  the  previous  analysis.  The  same 
methods  are  applicable  to  the  estimation  of  other  spectral  parameters. 

2 

The  estimates  of  f^  and  h  can  be  realized  by  determining  the  coefficients  c^  or  a 
spectral  estimate  <fr*(f)  by  some  standard  method  and  then  carrying  out  the  search  for 
the  maximum  likelihood  estimate  on  a  digital  computer.  An  explicit  interpolation 
method  such  as  Swerling’s  appears  practical  only  when  the  parameters  are  known  to 
lie  within  a  range  which  is  narrower  than  the  width  of  the  principal  fluctuations  of  the 
spectrum. 
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